Theory and Simulation Framework for the Relaxation of Nuclear Spin Order in Porous Media

The theory of nuclear spin relaxation in a liquid permeating a solid structure of irregular geometry is examined. The effects of restricted diffusion and the demagnetizing field generated by an inhomogeneous distribution of magnetic susceptibility in the system are explored. A framework comprising Brownian Dynamics, average Hamiltonian theory, and Liouville-space spin dynamics is proposed for simulating nuclear spin relaxation in 3D models of random structures obtained from CT scans of actual samples. Simulations results are compared with experimental data. An analytical solution valid within approximation is also reported.


■ INTRODUCTION
Nuclear spins diffusing in solution experience a variety of relaxation mechanisms that limit the lifetime of nonequilibrium spin order. 1 Even the most fundamental interaction with a magnetic field can induce spin relaxation, when the field is inhomogeneous and the spins are free to diffuse in solution. 2,3 At the basis of this effect, there is a phenomenon well-known to scientists performing diffusion NMR or MRI experiments, i.e., diffusive attenuation. Molecular translational diffusion in a spatially inhomogeneous magnetic field leads to a phase shift of the transverse spin order. 4 This effect has been long exploited to measure diffusion coefficients and diffusion tensors. 5 Diffusive attenuation also limits resolution in magnetic resonance imaging and, more relevant to this work, can lead to an additional source of relaxation for transverse spin order. The situation is exacerbated when molecular diffusion takes place in the voids of porous media. The difference in magnetic susceptibility between the solid framework and the voids in such systems give rise to local inhomogeneties in the magnetic field. Once a porous medium is placed in a static magnetic field, and the molecules move within and across pores, the nuclear spins carried by such molecules experience a magnetic field that fluctuates in both magnitude and direction. The effect of a randomly fluctuating magnetic field on nuclear magnetization has first been described by Redfield 6 and Abragam; 1 Callaghan 4 provides a thorough review of the spin dynamics of molecules undergoing diffusion in the presence of magnetic field gradients.
For a spin diffusing in a porous system, the statistical properties of this random field depend both on the dynamics of diffusion and on the spatial inhomogeneity of the magnetic field, which are both influenced in a complex manner by the spatial structure of the pores. Quantitative prediction of the relaxation caused by this effect poses therefore several challenges: 1. measurement of the random spatial structure of the pores and representation by an appropriate statistical model; 2. calculation of the resulting inhomogeneous magnetic fields in the pore space; 3. modeling of the restricted diffusion 2 of molecules in the pores and computation of the fluctuating field as a function of time; 4. conversion of the fluctuating field's statistical properties into a quantum mechanical propagator that describes spin relaxation.
Various aspects of the effect of susceptibility inhomogeneities on NMR in porous media have been discussed in the literature, 7−15 with most studies driven by the importance of diffusion NMR and MRI in materials science, in the petroleum industry, or in medicine. The phenomenon has become particularly relevant for us in relation to our interests in long-lived spin order. 16 In recent studies, we have exploited the long position-encoding time offered by long-lived spin order in diffusion-NMR applications to localize spins over long times or long distances in magnetic resonance imaging, 17 measure very slow flow, 18 extend the scope of diffusive-diffraction q-space imaging, 19 or measure tortuosity in porous media. 20 The long-lived spin order is created from single quantum coherences, which are converted by special pulse sequences. While the long-lived order is unaffected by diffusive attenuation, this does not apply to the single quantum coherences it originates from. Creation and observation of long-lived spin order is severely limited in porous media at high magnetic fields. The pulse sequences, for example M2S/S2M 21 or SLIC, 22 must be synchronized to the actual spin system's parameters and therefore cannot be kept conveniently short to minimize the effects of diffusive attenuation. In particular, the M2S/S2M sequence relies on relatively long and repeated echoes, while SLIC is based on a long spin lock. In either case, coherence decay typically results in a less efficient interconversion between long-lived spin order and single quantum coherence. As a result, when these pulse sequences are run in certain porous media no conversion of transverse into long-lived spin order is obtained at all. This is somewhat bizarre considering the fact that the long-lived spin order itself is unaffected by field gradients in the sample (unless the field varies on a molecular length scale). Motivated by the necessity to characterize this phenomenon, we have developed a simulation framework that predicts the relaxation of spin states due to susceptibility inhomogeneities in a porous system of arbitrary complexity. It relies on microcomputed X-ray tomography scans (μCT) of the porous structure, from which the positiondependent field is computed by a discrete Fourier transform approach. The same structure is used for a Brownian dynamics simulation of diffusion in the liquid phase. The resulting trajectories then serve as input to an average Hamiltonian approach to compute the relaxation propagator. As we discuss in the following, diffusive attenuation can also be treated theoretically, but fairly strong assumptions on the stochastic independence of diffusion and the random magnetic fields must be made.
We tested the predictions from both simulation and theory against experimental data taken at three different magnetic fields on three model porous media made up by randomly packed polyethylene spheres of different size imbibed with a tetramethylsilane solution in methanol-d 4 .

■ THEORY AND SIMULATION METHODOLOGY
Magnetic Fields in Porous Media. The magnetic field at each point in the structure of a porous medium can be thought of as the superposition of the magnetic field generated by the dipole magnetic moments present at each point in the structure. In general, the magnetic induction B (which governs the quantum dynamics of nuclear spins) is related to the magnetic field H and the magnetization M as In diamagnetic and paramagnetic materials, the magnetization is proportional to the magnetic field H 0 as in In the present context, χ is a piecewise-constant function of position. If we assume a porous solid described by the structure function S(r), the function χ(r) is given by where χ s and χ L are the volume susceptibilities of the solid and liquid phases, respectively. Obviously, this continuum treatment is only applicable for structures that are much larger than atomic sizes. This is adequate in the present context, as we are dealing with micrometer-sized pores. Without loss of generality, we assign the z-axis to be aligned with the external magnetic field, such that H 0 = H 0 e z . The magnetization then becomes The magnetic field has to satisfy Maxwell's equations By inserting (1) into (5), we obtain The magnetic field can be separated into a constant (external) field, which we assume to be aligned with the z-axis, and a position dependent field H d (r) as in For historical reasons, H d is known as the "demagnetizing field".
In the presence of piecewise-constant magnetic susceptibility, H d varies continuously as a function of position. It is these variations that cause diffusive attenuation. Introducing (8) into (7), we obtain With (4), and noting that ∇·= (∂ x , ∂ y ,∂ z )·, this becomes Since the structure function S(r) is constant everywhere except on its boundaries, this differential equation, together with Ampere's law (6), represents an inhomogeneous boundary value problem, which can be solved using the Green's function r r G r e r r e ( ) 1 4 The demagnetizing field is then given by the convolution 23,24 H S H r G r r r r Note that μG(r − r′) represents the magnetic field generated at position r by a point magnetic dipole μe z located at r′. The demagnetizing field is therefore the result of local magnetic dipoles that are induced in the paramagnetic or diamagnetic material by the external field H 0 .
It will prove convenient to formulate the problem in dimensionless form. In this way, the results obtained can be scaled to different pore sizes, diffusion coefficients, and magnetic susceptibility differences. We therefore introduce a dimensionless position vector l r = (13) where l is a length scale which can be chosen freely. It could, for example, represent the average pore size, or the voxel size of the CT data. The Green's function G(r) is homogeneous of degree −3; i.e., it has the property The convolution (12) can therefore be rewritten as 15) This means that the demagnetizing field is length-scale invariant. Using the Fourier convolution theorem, the demagnetizing field can be computed as where is the three-dimensional Fourier transform, 1 is its inverse, and s(ξ) = S(lξ) is the dimensionless (or scale-free) structure function. This formulation is convenient for numerical computations, since G( ) [ ] need only be computed once, and results derived for one structure function s(ξ) are valid at any length scale.
Fluctuating Field Experienced by Diffusing Spins. As the molecules diffuse in the pore structure, they randomly sample the demagnetizing field, such that each diffusing molecule experiences a field that fluctuates randomly as a function of time. The statistical properties of these fluctuations arise as a combination of the statistics of the Brownian motion of the molecules and those of the demagnetizing field. These two random processes are potentially highly correlated�after all, the structure S(r) both constrains the diffusion and gives rise to the demagnetizing field.
The Brownian motion of a diffusing molecule makes its position a random function of time r(t). If the diffusion is free, then this function represents a Wiener process with where Δr = r(t) − r(0), and D is the free diffusion coefficient. This continues to be valid for the restricted diffusion inside the pore structure S(r) for diffusion times that are short compared to a 2 /6D, where a is a measure of the pore size. For longer times, collisions with the pore walls cause deviations from free diffusion, and Δr(t) is no longer a Wiener process. The motion of the particle inside the pore structure makes the magnetic induction experienced by the particle a stochastic function of time: Averaged Propagator. As a model system, we consider an ensemble of molecules containing a single nuclear spin with spin quantum number 1 / 2 . In addition to the fluctuating demagnetizing field B d (t), the molecules are exposed to an external magnetic field assumed along the z-direction of the laboratory frame, H 0 e z . Since the transverse components of the demagnetizing field give rise to terms in the spin Hamiltonian that do not commute with the Zeeman term arising from H 0 e z , an average Hamiltonian approach is required. 25 Average Hamiltonian. The total time-dependent Hamiltonian Ĥ(t) for a spin-1 / 2 system diffusing through a demagnetizing field and in the presence of a large external field H 0 e z is given by where 26 As mentioned above, the two terms in eq 19 do not commute with each other. However, in many cases of interest, the fluctuations of B d are on a much slower time scale than the period of the Larmor precession. Formally, this can be expressed by the condition In this case, it makes sense to work with the average Hamiltonian over a full Larmor precession. To do so, we transform the Hamiltonian into the interaction frame of the Zeeman term of the Hamiltonian, Ĥz, as in where ω 0 = −γ B 0 is the Larmor frequency. Equation 22 can be expanded by explicitly solving the frame rotation to obtain with ω Assuming the time dependence of the fluctuating field to be slow compared to the Larmor period, the time-dependent Hamiltonian above can be averaged over a Larmor period to give rise to the Magnus series: 27 where the terms on the right-hand-side are the first, second, third, and so on orders of the expansion. The first two orders are given by 25 And, for the current case we obtain In the context of this paper, only the first-order averaged Hamiltonian will be considered. This is because in diamagnetic and paramagnetic systems |χ s − χ L | ≪ 1.
Since ω x , ω y , and ω z scale linearly with the susceptibility difference χ s − χ L , the second order term is proportional to (χ s − χ L ) 2 and can therefore be neglected. Single Spin−Echo Propagator. The central task of our simulation framework is to numerically compute the averaged propagator for a spin-1 / 2 system, subjected to the first-order truncated averaged Hamiltonian in eq 27, during a single spin− echo pulse sequence consisting of a π x pulse separated by equal time delays τ (see Figure 1). The propagator for such a sequence is given by where Ĥ̂is the commutator superoperator of the Hamiltonian and T̂̂is the Dyson time-ordering superoperator. The first-order truncated average Hamiltonian Ĥ̂( 1) (t) commutes with itself at all times. Introducing the following notation for rotation superoperators: the propagator in eq 29 can be expanded to first order as where Û̂(0, 2τ) is the part of the echo propagator associated with the decay of transverse magnetization during the time interval 2τ. Using the explicit form of the average Hamiltonian in eq 27 and recalling the definition ω z (t) = −γB d,z (r(t)), we obtain with B d,z = μ 0 H d,z as calculated in eq 16. It can already be appreciated that the propagator Û̂(0, 2τ) will approach unity in the limit of echo times that are short compared to the correlation time of the fluctuating Larmor frequency ω z (t). The effects of the random field experienced by the diffusing spins can therefore be arbitrarily scaled down by applying π pulses separated by short enough time intervals. In practice, however, the available (and tolerable) RF power impose limits to how fast these π pulses can be repeated. Moreover, in several situations, such as in the case of pulse sequences that prepare singlet spin order, the echo time cannot be chosen arbitrarily short but it is rather dictated by the spin system parameters.
In an ensemble of molecules containing nuclear spins in a fluid phase permeating a porous solid, ω z (t) becomes a random variable, due to the diffusion of the particle within the pores. Hence, in order to describe the global evolution of the system, we need the ensemble average of eq 32, formally given by Relaxation of Longitudinal and Transverse Spin Order. The averaged propagator in eq 33 can be used to study the dynamics of the density operator during the pulse sequence in Figure 1. Still in the case of a single spin-1 / 2 system, assuming thermal equilibrium at the beginning of the experiment, the density operator immediately after the first 90 deg pulse is containing the transverse spin order, i.e., single quantum coherences. If we adopt the following basis set for the Liouville space: {Î−, Iẑ, 1, Î+}, then the density operator at time t = 0 is represented by the following column vector: The density operator at time 2τ is then calculated from In the same basis, the propagator Û̂e(0, 2τ) is represented by the matrix and, R̂̂x(π) is represented by Hence, we find Giving the fact that in quadrature detection we only observe the −1Q coherence we can restrict our attention to the term The ensemble averaged exponential coefficient describes the decay of the transverse magnetization as the spins diffuse in the inhomogeneous magnetic field. The factor can be numerically calculated as explained below. Incidentally, it is easy to see that this phenomenon applies to transverse but not to longitudinal spin order, and therefore, it contributes to T 2 but not to T 1 relaxation. In the basis adopted, the longitudinal spin order is represented by the column vector which is clearly left unaltered by the averaged propagator in eq 37 and simply switched in sign by the propagator in eq 39. Indeed, a signature of the presence of susceptibility inhomogeneities induced relaxation is the experimental observation that T 1 remains unaffected, unlike T 2 . Note that this fact is not dependent on the length scale of the susceptibility inhomogeneities, but follows purely from the average Hamiltonian (27) and (28). 2 Analytical Solution. Let us consider the propagator element P 11 in eq 38, which governs the evolution of observable transverse spin order. For small enough diffusion displacements one can use a first-order Taylor expansion. It is convenient to choose the middle of the spin echo interval (the time of the inversion pulse) as the expansion point and the time origin t = 0. This gives where ∇B z is the field gradient experienced by the diffusing spin at time t = 0. This leads to Since the diffusion displacements Δr before and after t = 0 are statistically independent, this may be simplified to The integral in eq 47 is the average diffusion displacement: With this substitution, we obtain To compute the ensemble average, we write the joint distribution function of the magnetic field gradient and the average diffusion displacement as The ensemble average of any quantity A can then be expressed as the integral As outlined before, the distribution function Φ(∇B z , R; τ) is unknown, and in general, it contains complex information on the interplay between the fluctuating field and restricted diffusion. A numerical approach to sample the distribution function Φ is described below. An analytical solution requires some simplifying assumptions, which will have to be justified in view of experimental results. If we assume the average diffusion displacement and the magnetic field gradients to be uncorrelated, the distribution function factorizes into This is a rather drastic assumption, which is only valid if the diffusive attenuation is so efficient that relaxation is essentially complete on the time scale that is needed for the molecules to reach the walls of the pores. We would therefore expect this assumption to be valid at high magnetic fields, where the magnetic field variations are large, and for large pore sizes. R is given by the integral (48). If we further assume that over short enough times τ the diffusion remains unrestricted by the pore walls, Δr is a three-dimensional Wiener process, with components where W t is the standard Wiener process with variance σ 2 (W t ) = t, and likewise for the other two components. The integral (48) can therefore be evaluated as where N(0,1) represents a Gaussian random process with zero mean and unity variance. We therefore have The matrix element P 11 can therefore be expressed as The spatial integral is in fact the three-dimensional Fourier transform of a Gaussian function, and can be evaluated analytically. This leads to This expression can be used to compute the propagator if the field gradient ∇B z is known. In the special case of ∇B z = G = const., this reduces to The Journal of Physical Chemistry B pubs.acs.org/JPCB Article which is the well-known echo attenuation coefficient derived by Carr and Purcell for a particle diffusing in a constant magnetic field gradient. 4,28 The ensemble average can be evaluated numerically as a volume average over the liquid phase. Using our earlier definition of the structure function S(r), which is unity in the solid and vanishes in the liquid, we obtain where ϕ L is the volume fraction of the liquid. The exponential in eq 57 can be expanded in a Taylor series. To leading order, this yields The field gradient scales proportionally to the susceptibility difference. Since the demagnetizing field is independent of length scale for a given spatial structure, the field gradient is inversely proportional to the pore size a. We can therefore set which remains valid under the assumptions discussed above, as long as the right-hand side is much smaller than unity. The characteristic attenuation time is given by

■ NUMERICAL SIMULATION FRAMEWORK
To numerically evaluate the propagator element P 11 in eq 38, we have set up a simulation framework (available in Mathematica 29 and Julia 30 languages) which consists of the following steps: 1. Digital rendering of the porous structure. A microcomputed tomography 3D image of the porous medium is acquired and binarized to obtain the structure function, s(ξ). This function contains information about the structural nature of the porous structure at each position ξ, namely s(ξ) = 1 if ξ falls within the solid matrix and s(ξ) = 0 if ξ falls within a pore. Note that such knowledge of the structure is discrete and has the same resolution of the CT scan (see Figure 2a) . The digitized CT structure is then Fourier transformed to yield s( ) [ ] which is required in eq 16 to calculate the field H d (r).

2.
Compute the magnetic f ield within the porous medium. The demagnetizing field H d (r) is computed on a 512 3 mesh of the CT image (taken at the center of the original image) using a discrete version of eq 16. The demagnetizing field is then linearly interpolated to create a finer resolution. 3. Simulate Brownian dif fusion. Using Brownian dynamics, the discrete translational trajectory resulting from the molecular random walk within and across the pores of the structure is calculated as follows: (a) An initial point ξ 0 = {x 0 , y 0 , z 0 } is randomly chosen to fall within a pore in the structure, i.e., s(ξ 0 ) = 0 and a step counter j is set to 1.
(b) A random step is attempted by setting ξ j = {x j−1 + δx j , y j−1 + δy j , z j−1 + δz j }, with the step increments δx j , δy j andδz j chosen randomly from a Gaussian distribution of zero mean and standard deviation Dt 2 s = , where D indicates the isotropic selfdiffusion coefficient of the molecule carrying the spin and t s is the time increment.
(c) The point ξ j is checked such that if s(ξ j ) = 0, then the new point is still in the pore space, the position ξ j is stored in the trajectory array p = {ξ j }, and j is incremented by 1. If s(ξ j ) = 1, a new random step is attempted by returning to step b.
(d) Steps b and c are repeated N times. This produces a discrete trajectory of N s entries (in general, N s ≤ N, due to the rejection of steps colliding with the solid phase). To produce a statistically accurate description of spins randomly diffusing in porous media, the procedure above must be repeated for a (large) number of initial positions so as to simulate the different translational trajectories accessible to the different molecules. The total time for each trajectory N s t s must be on the order of the experimental time of interest, typically in the range between tens and hundreds of milliseconds in a liquid-state NMR experiment. In this paper, we have adopted the alternative approach consisting of simulating a very long trajectory (hundreds of seconds) for a single molecule and dividing this into many shorter subtrajectories. The two approaches produce statistically equivalent results if the pore structure is sufficiently interconnected as we believe it is the case of our test samples below. The result of a typical random walk simulation  Table 1.
The Journal of Physical Chemistry B pubs.acs.org/JPCB Article performed on the test samples described below is shown in Figure 2c.

Compute f luctuating f ield experienced by dif f using spins.
Once the field map is generated (step 2) and molecular positions as a function of the time are computed (step 3), the fluctuating magnetic field experienced by the spins along their translational trajectory becomes available. The field experienced by the spin at the jth step of its random diffusion through the porous medium is Since each point in the trajectory is spaced in time by t s , the field experienced by the spin at the time t = n s t s (i.e., after n s steps) is given by where only the real part of the complex exponential function appearing in eq 38 needs to be calculated since the imaginary part vanishes due to the detailed balance. All parameters used in the simulations are collected in Table 1.

■ EXPERIMENTAL MATERIALS AND METHODS
In order to compare numerical predictions with experiments so to validate our simulation framework we have measured the signal area as a function of the echo time τ using the pulse sequence in Figure 1 with n = 1 and in a number of model porous media made up by polyethylene beads of different size. These measurements have been done at three different values of the static magnetic field. For completeness, we also measured T 1 and T 2 decay times in all samples described below. Samples. The model porous media samples used in this work were made by random packs of polyethylene beads. All samples were prepared in a 5 mm NMR tube sourced from Wilmad LabGlass. The total packing height was 5 cm in order to fully encompass the probe coil region with sufficient excess to ensure that the packing was as uniform as possible across the region of interest. The packing was done by weighing out ca. 1.2 g (bead size dependent) of the polyethylene beads and adding this to the NMR tube in between 2 and 4 aliquots, depending on the bead size. Between each aliquot gentle manual tapping was undertaken to aid in packing. Polyethylene microspheres were purchased from Cospheric's CMPS products in sizes 212−250, 500−600, 710−850, and 1000−1180 μm and were used without further purification. The bead packings were imbibed with a 0.193 M stock solution of tetramethylsilane (TMS) in deuterated methanol (MeOD). Both solvent and solute were purchased from Sigma-Aldrich and used without further purification. Where required, gentle manual shaking was utilized to remove any visible air bubbles. Table 2 summarizes the sample composition and nomenclature.
NMR. NMR data were collected at three different values of the static magnetic field, namely 7.05, 9.4, and 16.45 T. Data at 7.05 T was collected in a Bruker Avance III 300 MHz spectrometer running topspin 3.5 and equipped with a 10 mm MICWB40 Bruker probe. Data at 9.4 T was collected on a Bruker Avance Neo 400 MHz spectrometer running topspin 4.0.8 and equipped with a 10 mm BBO Bruker probe. Data at 16.45 T was collected on a Bruker Avance Neo 700 MHz spectrometer running topspin 4.0.7 and equipped with a Bruker CPP TCI 700S3 probe. For all data collected, a single shim file was saved based on shimming the beads-free BLK sample (see Table 2). The same sample was also used to optimize the 90 and 180 deg pulses at all fields. Data fitting and relaxation values were calculated using a custom-made palette running within the Wolfram Mathematica software. In the following, we will refer to T 1 , T 2 , and single echo experiments. T 1 measurements use a standard saturation recovery pulse sequence. T 2 experiments use a standard CPMG pulse sequence like the one reported in Figure  1 run with fixed τ and variable n. Single echo experiments are also based on the pulse sequence in Figure 1 but are run with n = 1 and at variable τ values.
μCT. To obtain the digitized structure function S(r), we used microcomputed tomography (μCT). CT images of the three model porous media in Table 2 were collected using a modified 225 kVp Nikon/Xtek HMX scanner. To improve contrast, the sample packings were prepared as described above but in   The Journal of Physical Chemistry B pubs.acs.org/JPCB Article custom-made 10 mm poly(carbonate) tubes formed from a cut section of 10 mm outer diameter and 1.5 mm wall-thickness tubing (Clear Plastic Supplies, U.K.). For all samples, 3 g of microspheres was added to the tube in a number of small aliquots. At each addition, the sample was gently tapped to improve the packing. The packed beads were imbibed with a 0.193 M solution of TMS in MeOD and manually shaken to remove any visible air bubbles. The tube was sealed and left overnight before imaging. μCT images were processed using the ImageJ software. Raw data comprising the whole sample where cut to manageable sizes for simulations (512 × 512 × 512 pixels taken from the center of the reconstructed volume). These reduced-volume images were filtered and binarized. A slice of each data set is shown in Figure 3. Binarization was achieved using ImageJ's automatic threshold tool with thresholds calculated for each slice and manually checked after binary conversion. A median filter of 3 pixels was applied to the volume before binarization in order to minimize noise and remove artifacts. A further median filter of 2 pixels was applied following binarization.
■ RESULTS AND DISCUSSION 1D NMR Spectra. 1D NMR spectra of the four samples in Table 2 have been recorded at 16.45 T to examine the effect of the susceptibility inhomogeneities on the spectral line width. The spectra are compared in Figure 4. The peak at 0 ppm is due to TMS while the peak at 4.8 ppm is due to the residual protonated methanol. As expected, the presence of the PE beads in the sample severely compromises the spectral resolution. The effect is more severe for the sample with small beads (PES) and lessens as the bead size increases. Clearly, the line width is dominated by the field inhomogeneity due to the susceptibility differences between the beads and the solvent. However, the magnetic field is scale-invariant, as discussed above. The fact that the line width does depend on the length scale, if only weakly, already suggests that there must be a contribution to the line width from a T 2 relaxation mechanism that becomes more effective at smaller length scales.
Relaxation Measurements. The T 1 relaxation decay constants of all samples were measured, at each magnetic field, using a conventional saturation recovery pulse sequence with the relaxation variable delay time ranging from 0.1 ms to 32 s. The results of these experiments are reported in Table 3.
As expected, the beads and resulting susceptibility inhomogeneities do not alter T 1 . Conversely, as established in the theory section above, the relaxation decay constant of transverse spin order, T 2 , measured using a train of spin echoes (Figure 1 with fixed τ and variable n�a.k.a. the Carr−Purcell−Meiboom−Gill method) must depend on the actual echo time chosen in the experiment. To highlight this effect, we have measured the T 2 decay constants in the PEM sample, at the three different magnetic fields and for a series of echo times, as indicated in Table 4.
The data in Table 4 show how the measured T 2 value gets shorter and shorter as the echo time is increased. This is because the spins can diffuse further away from their original positions during an increased echo time. The values of T 2 also become consistently smaller as the static magnetic field is increased as the effect is directly proportional to the magnetic field (see eq 41).
Diffusion Experiments. The value of the unrestricted diffusion coefficient of TMS in MeOD is required in the random walk part of the simulation framework. This coefficient has been experimentally measured at 7.05 T using sample BLK (see Table  2) and a conventional pulsed gradient stimulated echo (PGSTE) sequence that used bipolar gradients and one spoiler   Table 2). Spectra b, c, and d are normalized to the same intensity scale for comparison.  Single-Echo Experiments. The results of single-echo experiments, using the pulse sequence shown in Figure 1 with n = 1, are compiled in Figure 5a for different magnetic fields and Figure 5b for different bead sizes. All decay curves exhibit a similar shape, reminiscent of a Gaussian function. At short echo times, the normalized signal area is constant, and starts to decay only once the echo time reaches several milliseconds. Clearly, smaller beads and larger magnetic fields lead to faster decay.
The scaling law, eq 61, predicts a decay with 1 − P ∼ τ 3 . Figure  6 shows a doubly logarithmic plot of the experimental results. In close agreement with the theoretical prediction, the data present as an initial straight line with a slope close to 3. At longer echo times, the decay slows down with P approximating zero.
Agreement with the theory developed above is further corroborated by plotting all available data as a function of the reduced echo time att . The characteristic attenuation time has been calculated according to eq 62, with the bead diameter used for the characteristic length scale a. This is shown in Figure 7. All data points, taken with different bead sizes and different magnetic fields, collapse onto a single master curve. This behavior is directly predicted by eq 57. The experimental data therefore directly validate the assumptions made in developing the theory, first and foremost the stochastic independence of the diffusion process and the irregular magnetic field distribution.
However, the theory only provides the scaling behavior of the propagator, not its absolute values, since the length scale a is not uniquely defined. For this reason, numerical simulations have been performed on all model porous media samples in Table 2 and at each magnetic field for which experiments were available. The agreement between simulations and experiments is quite good considering that we do not know the exact magnetic susceptibilities of the samples. The deviation from experimental measurements increases at large value of τ. This can be related to a gradually more important contribution of other relaxation mechanisms that appreciably contribute to T 2 at larger τ. It is also noteworthy that the diffusive attenuation effects are less severe as the beads sizes (and therefore the pore size) is increased. This is due to the fact that for large pore size the spin have to diffuse much further before to appreciate the fluctuation nature of the demagnetizing magnetic field due to susceptibility inhomogeneities. In parts d−f of Figure 8, we have compared simulations and experiments for the model porous media sample containing the smallest beads (and therefore the smallest pore sizes) for three values of the magnetic field, namely 7.05, 9.4, and 16.45 T, from parts d to f, respectively.
Once again, we note a very good agreement between simulations and experiments which confirm the good quality of our simulation approach. As expected, the diffusive attenuation effects are more severe as the field is increased because the demagnetizing field depends linearly on the static magnetic field.
In order to estimate at which magnetic field these effect becomes negligible at long echo times that would be required, for example, in a M2S/S2M experiment, we have simulated the effects of diffusive attenuation for the PES sample at 3 smaller values of the static magnetic field, namely 1, 0.5, and 0.1 T. The results of these simulations are reported in Figure 9.
It becomes clear that the diffusive attenuation at echo time duration of, say 30 ms, requires the experiments to be taken at magnetic fields of 0.1 T or lower when molecules are diffusing between the pores resulting from the packing of 212−250 μm plastic spheres, for which a susceptibility mismatch of 2.71 ppm exists between the structure and the solvent. The simulation framework is therefore very useful in estimating these properties for any actual porous structure providing that the involved magnetic susceptibilities are known and an accurate image of the sample can be acquired either through CT scan, as in this work, or via MRI and indeed any other available techniques.  Table 2). The signal area has been plotted so that the shortest time point has intensity 1.   Table 2). The signal area has been normalized so that the shortest time point has unitary area.

■ CONCLUSIONS
A theoretical framework has been developed to predict the diffusive attenuation of transverse nuclear spin order in liquids imbibed into random media. Experimental measurements have been performed on densely packed polyethylene beads of varying size and at a variety of magnetic fields. The scaling of the resulting attenuation with bead size as well as with the magnetic field agrees closely with theoretical predictions derived by neglecting correlations between the susceptibility inhomogeneities and molecular diffusion. In order to provide a quantitatively accurate prediction, a simulation framework has been developed that departs from CT images of the 3D structure of the sample. The CT data is used to compute the inhomogeneous demagnetizing field, and as a framework for Brownian dynamics simulation of molecular diffusion. This provides time-series data for the magnetic field experienced by a diffusing molecule, from which the required ensemble averages can be computed. The predictions were found to be in excellent agreement with experimental observations. The theoretical and simulation framework presented here enables systematic predictions of the diffusive attenuation behavior of more complex spin states. In turn, this is of great importance in the design of experimental strategies for magnetic resonance techniques based on long-lived spin states and their applications to study random media. As the simulation framework does not neglect correlations between magnetic field and diffusion, it may also be applied to low and intermediate field regimes where such effects are expected to be relevant. Actually, such experiments are currently underway in our laboratory and will be reported separately.  Table 2). The signal area has been normalized so that the shortest time point has unitary area. Simulations use the parameters in Table 1 and are the result of eq 65.  Table 1 and are the result of eq 65.
The Journal of Physical Chemistry B pubs.acs.org/JPCB Article